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ABSTRACT 

In this paper, we use large P 3 M N-body simulations to study the three-point cor- 
relation function £ of clusters in two theoretical models. The first model (LCDM) is 
a low-density flat model of Slo = 0.3, Ao = 0.7 and h = 0.75, and the second model 
(PIM) is an Einstein-de-Sitter model with its linear power spectrum obtained from 
Q\ observations. We found that the scaled function Q(r, u, v), which is defined as the ratio 

of C(r, ru, ru + rv) to the hierarchical sum £(r)£(ru) + £(ru)£(ru + rv) + £(ru + rv)^(r) 
(where £ is the two-point correlation function of clusters), depends weakly on r and u, 
but very strongly on v. Q(r, u, v) is about 0.2 at v = 0.1 and 1.8 at v = 0.9. A model 
^"~5 of Q(r, u, v) = 01O 1 ' 3 " can fit the data of ( very nicely with k, 0.14. This model is 

qq found to be universal for the LCDM clusters and for the PIM clusters. Furthermore, 

£Sj Q{ r , u , v ) is found to be insensitive to the cluster richness. We have compared our N- 

body results with simple analytical theories of cluster formation, like the peak theories 
or the local maxima theories. We found that these theories do not provide an adequate 
description for the three-point function of clusters. We have also examined the obser- 
vational data of £ presently available, and have not found any contradiction between 
\Q the observations and our model predictions. The independence of q in a projected cat- 

(~*> alogue of clusters is shown to be much weaker than the ^-dependence of Q found in the 

three-dimensional case. This is probably the reason why the ^-dependence of Q has 
not been found in an angular correlation function analysis of the Abell catalogue. The 
^-dependence found in this paper might be an important feature of clusters formed 
Q | in the Gaussian gravitational instability theories. Therefore it would be important to 

I search for the ^-dependence of Q in redshift samples of rich clusters. 

O 

%—{ Key words: galaxies: clustering - large-scale structure of the universe - cosmology: 

C/2 observations 



1 INTRODUCTION 

The observational fact that the two-point correlation function of rich clusters of galaxies has an amplitude much higher than 
that of galaxies (Bahcall & Soneira 1983; Klypin & Kopylov 1983), has been explained by Kaiser (1984) as a consequence of 
the hypothesis that rich clusters form at high peaks of the underlying Gaussian fluctuation field on the cluster mass scale. In 
his paper, Kaiser derived a relation between the two-point correlation function £ p k of peaks and the two-point function f m 
of the underlying mass. This derivation has been refined and generalized to the three-point correlation function £ p k of peaks 
by many authors (e.g. Politzer & Wise 1984; Matarrese, Lucchin & Bonometto 1986; Jensen & Szalay 1986). For a Gaussian 
perturbation field and in the limits of high sharp peaks and weak underlying clustering (£ m <C 1), the three-point function 
C P k is found to be 

( P k(ri2, r 2 3, rzi) = Qh ^ P k(r 12 )^pk(r23) + ipk(r23)i P k(r3i) + £pfc(r 3 i)£p;t(ri2)j + Qkipk(ri2)ipk(r23)ipk(r3i) (1) 

with Qh = Qk = 1. For this case, £ p k is said to have the Kirkwood form. If Qk is zero, £ p k is said to have the Hierachical 
form. Szalay (1988) relaxed the limit of high sharp peaks and examined £ p k for a general threshold function. Under the limit 
of weak underlying clustering, he showed that the cubic term Qk£pk{ri2)£pk{'r23)£pk{'r3i) always accompanies the hierarchical 
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term Qh[^ P k(ri2)^ P k(r2s) + £ P k(r23)£ P k(r3i) + £ P k(r3i)£ P k(ri2)] with Q k = Q 3 h even for a general threshold function (though 
some other terms may appear in this case). He suggested that this relation can be used to test the hypothesis that clusters 
form at high peaks of the primordial density field. 

Studies based on the projected catalog of Abell clusters, however, indicate that the three-point correlation function 
of clusters approximately obeys the hierarchical model with Qk = and Qh ~ 0.6 (Jing & Zhang 1989; Toth, Hollosi & 
Szalay 1989). The same hierarchical model seems to hold for of clusters selected from the Lick catalog (Borgani, Jing 
& Plionis 1992). A few workers have tried to extract the three-point function of clusters from redshift samples of clusters, 
however the constraints on the form of from these studies are rather weak because the redshift samples are too small (Jing 
& Valdarnini 1991; Davies & Coles 1993). Recent analyses of the skewness Ss(R) for the Abell clusters and the APM clusters 
seem to support the hierarchical model because the measured Ss(R) depends very weakly on R (Plionis & Valdarnini 1995; 
Cappi & Maurogordato 1995; Gaztanaga, Croft & Dalton 1995 hereafter GCD95). 

In real observations of £ c ;, as Coles & Davies (1993, hereafter CD93) pointed out, the condition f m <C 1 can be hardly 
satisfied because the two- and three-point correlation functions of clusters are very difficult to measure on the scales of f m <C 1. 
They therefore, using numerical integrations, investigated the form of £ p k as a function of {pi; for the quasi-linear regime f m < 1. 
They find that the hierarchical model is a better fit to £ p k than the Kirkwood model though neither model provides a good fit 
to C P k. As they pointed out, one problem with this peak theory is that a region with density exceeding some threshold may be 
identified with several distinct clusters because of its finite spatial extent, thus overestimating the correlations on small scales 
(Coles 1989; Lumsden, Heavens & Peacock 1989). They therefore studied the three-point correlation of local maxima of the 
density field. It is found that the hierarchical form can approximately describe the three-point correlation function of local 
maxima. However, because of the technical complexity of this approach, their conclusion is drawn only for the one-dimensional 
density field. One common prediction of both theories is that the amplitude Q of the three-point function increases with the 
density threshold. 

The main weakness of these analytical theories is that they are unable to treat the non-linear density fluctuations, 
especially merging of clusters. As pointed out by Croft & Efstathiou (1993) based on N-body simulations, these theories 
predict too high a two-point correlation function for clusters. This weakness should also influence the calculations of the 
three-point correlation function. Therefore it is obviously very important to use N-body simulations to calculate the three- 
point correlation function for clusters. Such a study is not only useful for testing the analytical predictions but also provides 
a basis to compare observations with theories. 

Here we report our results of the three-point correlation function £ of clusters in N-body simulations. The simulations 
are a set of P 3 M simulations of 10 6 particles in a simulation box of 300 or 400 /t -1 Mpc. A description of these simulations 
will be given in Section 2, where identification procedures of clusters will be discussed. In Section 3, we will report the three- 
point correlation function of simulated clusters, with emphasis on the dependence of £ on sizes and shapes of triangles. We 
will find that neither the hierarchical nor the Kirkwood model can provide a reasonable fit to the three-point function of 
simulated clusters. There we suggest another model which can fit £ quite well. Then in Section 4, we compare our results with 
previous works: observational statistics, analytical theories, and simulation studies of the skewness *. Our main conclusions 
are summarized in Section 5. 

2 N-BODY SIMULATIONS 

The simulations used here are two sets of large P 3 M N-body simulations. The first set (hereafter LCDM) simulates a low- 
density flat universe of fio = 0.3, h = 0.75 and as = 1, where fio is the current density parameter, h is the Hubble constant in 
unit of lOOkms^Mpc" 1 and as is the linearly evolved rms mass fluctuation in a sphere of 8 ft 1 Mpc at the present time. For 
this simulation, we assume that the primordial power spectrum is a Harrison-Zel'dovich spectrum. Its linear transfer function 
is taken from Bardeen et al (1986). The second set (hereafter PIM) simulates an Einstein-de-Sitter universe with as = 0.8 and 
a phenomenological linear power spectrum which was obtained by Branchini, Guzzo & Valdarnini (1994) from an analysis 
of redshift surveys of galaxies. Because of these parameters, both models provide good approximations to the real universe, 
although, as Branchini et al. noted, the latter model slightly lacks large-scale clustering power when compared with the APM 
and EDSGC angular two-point functions of galaxies (Maddox et al. 1990; Collins, Nichol & Lumsden 1992). 

The standard P 3 M technique is used for these simulations (Hockney & Eastwood 1981; Efstathiou et al. 1985; Valdarnini 
& Borgani 1991; Jing & Fang 1994). For both models, 10 6 particles are used in each simulation and 8 realizations are run for 
each model. For the LCDM simulation, we use a box of L = 400/t -1 Mpc and the force resolution r\ = 0.2/t -1 Mpc; for the 
PIM simulation we use L = 300 /t _1 Mpc and r] = 0.3/t _1 Mpc. Since we have simulated a huge volume of the model universes 
( 27 or 64 xlO 6 h~ 3 Mpc 3 ) with a large number of particles, we are able to calculate accurately the three-point correlation 

* Because the skewness depends on the three-point correlation through an integral, the shape dependence we find in this paper cannot 
be found in the skewness analysis 
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function of clusters in these models. 

Clusters are selected according to the single criterion — the mass overdensity T in a spherical volume of radius r c = 
1.5/t -1 Mpc. Similar criteria have been used for cluster identification in real observations (Abell 1958; Dalton et al. 1992; 
Lumsden et al. 1992). In order to study the richness dependence of the cluster three-point correlation function in detail, 
we have selected both rich and poor clusters with T > T mtn . T mtn is taken to be 90 for the LCDM model and 38 for the 
PIM model. With these overdensity thresholds, we can identify ?S 3200 clusters in each realization of the simulation for both 
models. Moreover each cluster in both models contains at least 20 particles. The poorest clusters in our simulations may be 
slightly poorer than the clusters in the above mentioned observations. We have used two methods to identify these overdensity 
regions (i.e. clusters). The first method is based on the friends- of- friends algorithm and a description of this method can be 
found in Jing et al. (1993; see also White et al. 1987). The linking parameter b in this work is taken to be 0.1 times the mean 
particle separation. This method has been applied only to the LCDM simulation because it is too CPUtime-consuming for 
large simulations (10 6 particles here). Therefore we developed another method which can work much faster. 

In the second method, we place a grid of Ng cells in the simulation volume. N g is 256 for the LCDM and 192 for the PIM 
simulations, so that the cell size in both simulations is 1.56 /t -1 Mpc. We count the number jC of particles in each cell. If a cell 
has the count p larger than the count in every of its 26 neighboring cells, this cell is recognized as a local density maximum 
and the cell center is regarded as the position of the maximum. For each local maximum, we calculate the barycenter Tq and 
the count Co of the particles inside a sphere of radius r c = 1.5/t -1 Mpc centered at the maximum. Around the barycenter 
Tq, we draw a new sphere of r c and calculate the barycenter r\ and the count C\ for the particles in this sphere. We repeat 
the same calculation for each updated barycenter, until the new count C t is not larger than the previous one C % -\. The 
previous barycenter r\_ 1 and count C,-i are then the center and particle count of a new cluster if its overdensity is larger 
than the threshold T mtn . This iterative calculation is attempted to minimize the effect of the uniform grid cells. Some clusters 
identified in this way may overlap, i.e. the separation between two clusters may be smaller than 2r c . We correct this simply 
by eliminating the smaller one of two overlapping clusters. A sample of clusters has thus been constructed. The cluster sample 
may still slightly depend on how the grid is placed. Especially some clusters may not appear as local maxima in this particular 
grid placement, thus they are not included in the sample. In order to reduce this effect, we displace the grid by 1/2 cell size 
along the x-axis and/or the y-axis and/or the z-axis (i.e. 7 displacements), and produce, in the way just described, 7 samples 
of clusters using these displaced grids. Although most clusters (~ 90%) are common to all the 8 samples, a small fraction of 
clusters appear only in some of these samples. Therefore we merge the 8 samples, and eliminate the overlapping clusters in 
the same way as we do for a single sample. This merged sample of non-overlapping clusters is our final sample of clusters. 

We have applied the second method both to the LCDM and to the PIM simulations. In order to test whether cluster 
properties, especially the results of this work, depend on the method of cluster identification, we have compared the mass 
functions and the two-point correlation functions of the LCDM clusters identified by the two different methods. The two 
identification procedures produce a consistent result for both measures. In Section 3, we will further show that these procedures 
give a consistent three-point correlation function of clusters. In the rest of the paper, unless explicitly stated, clusters used 
for our analysis are those identified by the second method. 



3 THE THREE-POINT CORRELATION FUNCTION OF CLUSTERS 

In order to study the richness dependence of the three-point correlation function, we construct 4 subsamples of clusters for 
each model according to cluster richness. The i-th subsample (i = 1, 2, 3, 4) is made of the 400 x i most massive clusters in each 
realization of the simulation, thus contains a total of 3200 x i clusters. For simplicity, we will call the i-th subsample of the 
LCDM clusters the LCDMi sample and the i-th subsample of the PIM clusters the PIMi sample. The overdensity thresholds 
r m8n of the LCDM samples are respectively 90, 136, 195, and 260, and the T mtn of the PIM samples are respectively 38, 59, 
90, and 125. 



3.1 The method of estimating the three-point correlation function 

The three-point correlation function C( r i2, ^23, J"3i) is defined by writing the joint probability dP\ 23 of finding three clusters 
concentrated in each of three volume elements dV\ , dV2 and dV 3 separated by r\ 2 , r 23 and r 3 \ as, 

dP 123 = n 3 [l + f (n 2 ) + f (r 23 ) + f (r 3 i) + ((r 12 , r 23 , r 3 i)] dVt dV 2 dV 3 (2) 

where n is the mean number density of clusters. Usually it is convenient to use variables r, u, v instead of r\ 2 , r 23 , r 3 \ (Peebles 
1980). For ri2 < r 23 < r 3 i, the relations of these variables are given as following: 

r 23 r 31 — r 23 . . 

r = r\ 2 , u = , v = . (3) 

ri 2 r 12 
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Figure 1 — The two-point correlation functions of four samples. The solid curves are the best fits of the models described in the text. 



The count N T (i,j, k) of cluster triplets with sides in the k)-bm R, < r < R,+i, Uj < u < + i and Vk < v < Vk+i are 

N T (i,j,k) =N (1) (i, h k) + N (2) (i, h k) + N (3) (i, h k) 

r R i+i r u j+i r v k+i 
N- ' (i, j, k) =8ir n N c i I I I r (u + v)u dr dn dv 

J R, Ju 3 JV k 

r R i+i r u ]+i r v k+i r 1 (4_'\ 

N {2 \i,j,k) =8ir 2 n 2 N c i ' ' ' , , . . 

N (3) (i, h k) =8Tr 2 n 2 N c i 



r (u + v)u \£,(r) + £,(ru) + ^(rti + rv) dr dn dv 



R, 



r 5 (u + v)u ((r, ru, ru + rv) dr du dv 



where N r j is the total number of clusters in the sample. For the N-body simulations where the periodic boundary is assumed, 
the above equation is strictly valid for the triplet count if the periodic effect is taken into account (in contrast with real 
observations where the boundary effect is more difficult to correct). If the two-point correlation function and the triplet count 
N T (i, j, k) have been measured, the three-point correlation function £(r, u, v) can be determined through Eq.(4). 

We therefore first measure the two-point correlation function £ for each sample. In Figure 1, we plot, as examples, the 
two-point correlation function for four samples. The error bars shown in the figure are the bootstrap errors estimated by 
the analytical formula of Mo, Jing & Borner (1992). For the LCDM clusters, we find that the two-point correlation function 
can be accurately fit by 

1 



(jj.1.5 _|_ ^ r 2 _|_ CJ ,3.5 

Because the PIM model has relatively less power than the LCDM model, the following formula 
1 



+ br 2 + cr 6 



(5) 



(6) 



can nicely fit the two-point correlation function of the PIM clusters. The fitting parameters a, b, and c are determined for each 
sample by a least-square fit to £ -1 (r) with the weight ^a^ 1 . The fitted £ are shown as solid curves in Figure 1. The figure (and 
a check with the fit results of other samples not plotted in the figure) shows that these fitting formulae can accurately express 
the two-point functions of the LCDM clusters for r < r max x 45 A -1 Mpc and of the PIM clusters for r < r max x 35 ft -1 Mpc. 
These formulae will be used to predict N (2) (i,],k). We have compared our cluster-cluster £ of the LCDM model with that of 
Croft & Efstathiou (1993) of a similar model, and found a good agreement between the two studies. Our results also confirm 
the weak richness-dependence of the £ amplitude originally found by Croft & Efstathiou (1993). 
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We then count the number N T (i, j, k) of cluster triplets in each sample. Here we use the following bins of r, u, and v. 
R l <r < R l+1 [\og(R t+1 / Ri) = 0.1, i = 1, • • -,i max , i?i = 3.1 A _1 Mpc]; 

U 3 <u< U J + 1 [log([/ J + 1 /^) = 1/6, j = 1, • • -, 6, U t = 1]; (7) 
V k < v < T4+i [T4 + i -V k = 0.2, k = 1, • • •, 5, V x = 0]. 

Since the two-point correlation function { is detected and fit well only at r < r max , we limit our analysis to triangles with 
their longest sides smaller than r max . Here we take i max = 12 (r max = 39/t -1 Mpc) for the LCDM clusters and i max = 11 
(jma X = 31 /t -1 Mpc) for the PIM clusters. Because of the constraint r?_,\ < r max , the integral upper limits of u and v for the 
N^'\i,j, k) now become min(?7j+i, r max /r) and min(T4+i, r max /r — n). 
We will use the scaled function 

nil \ ^( r ' u ' v ) 

h!(r, u, v) - ^ r ^ rM ) + £( ru )£( ru _|_ rv ) _|_ £( ru _|_ r „)^( r ) ^ 

to express the three-point correlation function. For the hierarchical model of Q(r, u, v) is a constant; for the Kirkwood model 
(Qh ~ Qk in Eq. 1), Q(r, u, v) decreases rapidly with r for {(r) > 1. Since it is difficult to show the full (3-D) dependence of 
Q(r, u, v) on r, u and v by simple means (e.g. plots), in the following we will separately show the dependence of Q(r, u, v) on 
each variable by taking an average of Q(r, u, v) over the other two variables. The least-square technique is used to find these 
averages Q(y) ( y = r,u, v). For example, Q(r) is found by minimizing 

2 S^f NT (hhk) - N (1) (i, h k) - N {2) (i, h k) - Q(i)N w (i,j,k)\ 2 

v n t (h h k) 

where a N T(i,j,k) is the bootstrap error of N T (i, j,k) estimated by the analytical formula of Mo et al. (1992), and 
/•-R.+i r u j+i r v k+i 

N (i, j , k) = 8ir 2 n 2 N c i / / / r 5 (u + v)u[£(r)£(ru) + £(ru)£(ru + rv) + £(ru + rv)£(r)] dr du dv . (10) 

J R, J u 3 J v k 

This least-square technique will also be used to test the hierarchical model and the Kirkwood model of £ and to estimate 
their parameters Q, Qh and Qk- For the Kirkwood model, we minimize 

2 W ^ T (»,J,fc) ~ N (1 \i, h k) - N (2 \i, h k) - Q h N (i \i, h k) - Q k N (5 \i, h k) \ 2 

X L>\ <T N T(i,j,k) ) ( ' 

with 



N (5) (i, j, k) = 8ir 2 n 2 N c i / / / r 5 (u + v)u£(r)£(ru)£(ru + rv) dr du dv . (12) 

J R, Ju 3 JV k 

For the hierarchical model, we just minimize the \ 2 °f Eq.(ll) with Qh replaced by Q and with Qu set to zero. 
3.2 The results 

In Figure 2, we show the results of Q(r), Q(u) and Q(v) of the LCDM clusters. The most noticeable point is the strong 
dependence of Q(v) on v. Q(v) increases with v, and its value is about 0.2 at v = 0.1 and about 1.8 at v = 0.9. The dependences 
of Q on r and u are much weaker: within ~ la error bar, a constant value ~ 0.5 is acceptable for Q(r) and Q(u). The decline 
of Q(r) at large r is due to the facts that v of the triangles with large r should be small because of the longest side r^i < r max 
and that Q(v) is an increasing function of v. For the range of cluster richness studied in this work, these dependences of Q 
do not depend on the richness, as shown by different rows of the figure. The uncertainty of Q of the LCDM1 sample is large 
because the number of clusters in this richest sample is relatively small. 

The results of the PIM clusters are shown in Figure 3. The Q(r), Q(u) and Q(v) of the PIM clusters are nearly the same 
as those of the LCDM clusters. This might indicate that these Q-functions do not sensitively depend on the density parameter 
fio or the cosmological constant Ao (these two parameters of the LCDM and the PIM models differ significantly). Since the 
difference of the linear power spectra between the two models is not large, it is difficult to say whether these functions depend 
on the power spectrum or not. 

We have also done the same analysis for the LCDM clusters identified by the friends- of- friends algorithm (Method I). 
The results are shown in Figure 4. The three Q-functions of this figure are nearly the same as those in Fig. 2. This means 
that our results are robust, not depending on the method used to identify clusters. 

The weak r-dependence and the strong ^-dependence of Q means that neither the Kirkwood model (with Qh = Qk) nor 
the hierarchical model is a good model for the three-point correlation function of clusters (for the range of scales analyzed 
here). We propose another empirical model in which Q(r, u, v) depends only on v through the relation 

Q(r,u,v) = 01O 13 " 2 . (13) 



6 Y.P. Jing, G. Borner and R. Valdarnini 



. i i i I i i i I i i i I i i i I i i i _l_ I i i i | i i i | i i i | _i_i i i | i i i | i i i | i i i | 



i - 



; i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 : : 



> 

u 

Z 

u 2 
o 

O" 1 



1 - 



2 h 
1 




LCDM1 



LCDM2 



I j 1 1 1 1 1 



1 



:: I I I I 



I I I I I I I I I I I I I I I I I I I 



LCDM3 



* j i 1 I i i 



I 



I I I 1 1 I I | I I I | I I I | I I I 



LCDM4 



i i . i i 



~i i i i i i i i i i i i i i i i 



1 1 1 1 1 1 1 1 1 1 



1 1 



1 1 1 1 1 1 1 1 1 1 



! I I I 1 1 



I I I I I I I I I I 



i i 



I I i 



I i i i I i i i I i i i I 




I I I | I I I | I I I | I I I | I I I ) 



i i I i i i I i i i I i i i I i i i 



0.4 0.6 0.8 1 1.2 0.2 0.4 0.6 0.2 0.4 0.6 0.8 1 



log 10 r (h _1 Mpc) l°g 10 u 



Figure 2 — The averaged Q of the LCDM clusters as a function of r (the left panels), u (the middle panels) or v (the right panels). Each 
row of three panels shows the results of one sample, with the sample name given in the right panel. The solid curves in the left panels 
are Q(v) = 0.15 10 1 ' 3 " 2 . 



As shown by the \ 2 tests, this model (hereafter Model III) can fit the simulation data of ( much better than the Kirkwood and 
the hierarchical models. Table 1 lists Xmzn °f these three model fits to ( of the LCDM clusters (cf. Eq. 11). The Xmzn °f Model 
III is much smaller than those of the other two models, which means that Model III can better fit the data of £ than the other 
two models. For 239 to 240 degrees of freedom, the Xmm °f the LCDM4 model is 71 for Model III, 561 for the hierarchical 
model, 671 for the Kirkwood model with Qh = Qk, and 544 for the Kirkwood model with two free parameters Qh and Qu, 
therefore the hierarchical and the Kirkwood models fail to describe the simulation data. Since the absolute values of Xmzn 
depend on the error model (here the bootstrap error) used in the least-square fits, we caution readers against over-interpreting 
the statistical meaning of these absolute values, though the bootstrap error was shown to be a good error model for such fits 
(see Mo et al. 1992 for a very detailed discussion). The best values of from these fits are given in Table 2. is about 0.15 
for the LCDM clusters and 0.13 for the PIM clusters. The solid curves in Figs. 2-4 are the Model III predictions with these 
values. It is clear that Model III can very nicely fit the simulation data of 
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Figure 3 — The averaged Q of the PIM clusters as a function of r (the left panels), u (the middle panels) or v (the right panels). Each 
row of three panels shows the results of one sample, with the sample name given in the right panel. The solid curves in the left panels 
are Q(v) = 0.13 10 
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4 DISCUSSION 

First let us compare our N-body results with the predictions of some analytical theories of cluster formation (CD93). The 
analytical theories explore the clustering properties either of regions with density fluctuation above some threshold (hereafter 
the peak theory) or of local linear density maxima above some threshold (hereafter the maxima theory). The peak theory 
predicts that Q slightly decreases with the increase of v (cf. Figs. 1-2 of CD93), contrary to our N-body results. The maxima 
theory can predict the three-point correlation function only for the one-dimensional case (i.e. v = 1), so it is unknown how Q 
changes with v in this theory. Both theories predict that Q increases with the peak (maxima) threshold, which is inconsistent 
with our result that Q does not depend on the richness of clusters. The failure of the analytical theories to match our N-body 
results implies that these theories can not be used to predict the three-point correlation function for clusters. It is not very 
surprising that the analytical theories are inconsistent with our N-body results, since these theories have also met problems 
to predict the two-point correlation function of clusters (Croft & Efstathiou 1993). 

We have noted two recent works which have analyzed the three-point correlation function and/or skewness Ss(R) for rich 
clusters in N-body simulations. Watanabe, Matsubara, & Suto (1994) analyzed the three-point correlation function and Ss(R) 
for simulation clusters. However, their cluster samples are too small, so that they are unable to address the dependences of Q 
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Figure 4 — The same as Fig. 2, but the clusters are identified by the friends- of- friends algorithm. 

Table 1 X 2 mln °f three model fits to the LCDM samples 

Hierarchical Kirkwood 1 Kirkwood 2 Model III 
LCDM1 51.0 54.4 50.4 37.5 
LCDM2 80.0 94.1 75.4 43.4 
LCDM3 182.2 223.5 178.4 52.3 
LCDM4 560.8 671.1 543.8 71.4 

The number of data points for these fits are 241. In Kirkwood 
1, Q h = Q k is assumed. In Kirkwood 2, both Qh and are 
free fit parameters. 

Table 2 The fit values of Q for the LCDM 
and PIM clusters 



Sample No. 
1 

2 
3 
4 



LCDM 
0.16 ± 0.03 
0.15 ± 0.02 
0.15 ± 0.01 



PIM 
0.11 ± 0.02 
0.13 ± 0.01 
0.136 ± 0.005 



0.140 ±0.004 0.131 ±0.003 
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Figure 5 — The skewness S3 (R) of an LCDM cluster sample measured with the count-in-cell (CIC) analysis (squares connected with solid 
lines), compared with the integral results of Equation (14) for the Hierarchical model (dotted line) and for Model III (thick solid line). 
The value at R = 10 /i -1 Mpc given by the CIC analysis is slightly low due to the finite size of the clusters. The numerical integration 
accuracy of Ss(R) is controlled to ~ 10%. Within this accuracy, the integral results given by both models are in good agreement with 
the result from the CIC analysis. 



on r, u and v, though their average value of Q is consistent with our results. 

GCD95 have analyzed only the skewness of clusters (but not the three-point function) in a set of CDM simulations. The 
normalized skewness S3 (R) of a spherical volume is related to the three-point correlation function through 



S 3 (R) 



C(R) 



((R) \ dr-i dr 2 dr?X(ri2, r 2 s, r 3 i) , (14) 

sphere R 



=-L / dn dr 2 S,{r 12 ) 



V 2 , 

•/sphere R 

where V = 4wR 3 /3. We have measured 53 (i?) for our LCDM cluster samples by the count-in-cell analysis (Peebles 1980) 
with a result in good agreement with that of an LCDM model of GCD94. Our measured 53 (i?) is also nicely consistent, as 
expected, with the 53 (i?) numerically calculated through Eq.(14) using Model III. This is shown in Figure 5. 53 (i?) is almost 
a constant ~ 1.7 for R < 25 A -1 Mpc (we limited our analysis to R < 25 A -1 Mpc because the three-point correlation function 
is measured only for r?_,\ < 40/t -1 Mpc. The constant S?_,(R) might extend much beyond R = 25/t -1 Mpc). 

A constant Ss(R) was usually regarded as a support for the hierarchical model in the literature. Actually the hierarchical 
model with Q PS 0.55 can reproduce, using Eq.(14), our measured Ss(R) (see Fig. 5), but the hierarchical model clearly fails 
to account for the direct count N T (i, j, k). It is not surprising that the different models for ( can lead to the same result 
for 53 (i?), since Ss(R) is an integral of £. Clearly a constant Ss(R) is not sufficient to argue for the hierarchical model. The 
count-in-cell analysis (or the moment method) has been widely used to study high-order correlation functions of extragalactic 
objects. Although these analyses have already yielded useful information for high-order correlations (e.g. Plionis & Valdarnini 
1995; Cappi & Maurogordato 1995; GCD95), the strong ^-dependence of Q found in this work will be definitely lost in the 
moment method. 

The average value of Q (i.e. averaged over r, u and v) of our model clusters is about 0.55, consistent with the observational 
results from the statistical analyses of Abell clusters, Lick clusters and APM clusters (Jing & Zhang 1989; Jing & Valdarnini 
1991; Borgani et al. 1992, Toth et al. 1989; Plionis & Valdarnini 1995; Cappi & Maurogordato 1995; GCD95). Of these 
observational works, only JZ89 have examined the richness-dependence and the ^-dependence of Q. They analyzed the three- 
point function for Abell clusters of richness R > 1 and R > 2 respectively, and found that Q does not depend on the richness 
of clusters. Their result is consistent with our results here. As for the ^-dependence, because the Abell catalogue is a projected 
sample, they measured the value of q in two bins of < v < 0.5 and 0.5 < v < 1 [where q(9, u, v) is defined in the same way as 
Q but with the angular two and three-point correlation functions, see JZ89 for details]. They found that q does not strongly 
depend on v, at least much more weakly than the ^-dependence of Q found here. The projection effect in a two-dimensional 
catalogue may have reduced the ^-dependence. We have examined this problem using the Limber equation (Peebles 1980). 
Assuming that the three-point correlation function obeys Model III with O = 0.15, the two-point function £ oc r -2 and the 
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Figure 6 — The u-dependence of q predicted by the Limber equation for a projected catalogue of rich clusters, if the cluster three-point 
correlation obeys the model of this paper. 

cluster sample is a perfect distance-limited sample of redshift z = 0.2, we calculate q(9, u, v) for different values of 9, u and 
v. We find that q depends very weakly on 9 and u, but does depend on v. The ^-dependence of q is shown in Figure 6. We 
see that q increases from q fa 0.5 at v = 0.1 to q fa 0.9 at v = 0.9. The ^-dependence of q in the two-dimensional case is 
much weaker than that of Q in the three-dimensional case. From Figure 6, we expect that the value of q in the two bins of 
< v < 0.5 and 0.5 < v < 1 is about 0.5 and 0.7 respectively. These two values are consistent with the result of JZ89 within 
the statistical uncertainty of their work (Ag fa 0.2). So it is not surprising that they could not discover the ^-dependence of 
q in their work. Overall the observational results are not contradictory to our result based on N-body simulations. 

It is worth emphasizing that the analytical theories compared above are simple theories of cluster formation. Their failure 
to explain our N-body results only means the limitation of these simple theories. A similar ^-dependence to that found 
here has been predicted for the mass three-point function £ m of the initially Gaussian fluctuation field by the second-order 
perturbation theory (Fry 1984). Although the relations between the cluster £ and the mass £ m are yet unknown, it seems that 
the non-Gaussian underlying density field (which is resulted from the dynamical evolution of a Gaussian density field) has an 
important impact on the form of the three-point function of clusters (see Fry & Gaztanaga 1993 for a general discussion of 
the biasing issue in the second-order perturbation theory). This could also be the reason why the peak theories or the local 
maxima theories fail to explain our N-body results, since these theories assume that the underlying density field is Gaussian. 
The relations between £ and £ m will be explored in a subsequent paper (Jing et al. in preparation). At this moment, it is 
reasonable to argue that the ^-dependence of Q probably is an important feature of clusters formed in gravitational instability 
theories of Gaussian fluctuation fields. As an example of non-Gaussian models, the toy model of a bubble universe in which 
galaxies are distributed on shells and clusters form at the "knots" of three crossing shells predict a very weak ^-dependence of 
Q for rich clusters (Jing 1991). This toy model is inspired by a non-gravitational instability theory of galaxy formation — the 
explosion scenario (Ostriker & Cowie 1981). Therefore, it is worthwhile to conduct a statistical analysis for the ^-dependence 
of Q in redshift samples of rich clusters. The ^-dependence of Q could serve as a discriminator between different scenarios of 
cluster formation. 

5 CONCLUSIONS 

In this paper we have used two sets of N-body simulations to study the three-point correlation function of clusters in theoretical 
models. We found: 

(1) The scaled three-point function Q(r, u, v) depends weakly on r and v, but very strongly on v. This function is found 
to be universal for the LCDM clusters and for the PIM clusters. Furthermore, this function does not depend on the cluster 
richness. 
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(2) Two well-known models of £, the Kirkwood and the Hierarchical model, cannot be accommodated to these features of 
Q(r, u, v). We proposed another model Q(r, u,v) = Q 10 1 ' 3 " which can fit the data of £ very nicely with PS 0.14. 

(3) Simple analytical theories of cluster formation, like the peak theories or the local maxima theories, fail to explain the 
weak dependence of Q(r, u, v) on the cluster richness and/or the ^-dependence of Q(r, u, v) found in our N-body simulations. 
The reason might be that these theories cannot describe the non-linear density fluctuation, especially the merging processes 
of clusters and that these theories assume a Gaussian fluctuation for the underlying density field. 

(4) The ^-dependence of q in a projected catalogue of clusters is found to be much weaker than the ^-dependence of Q. This 
is probably the reason why JZ89 have not found the ^-dependence in their analysis of £ for the Abell catalogue. The weak 
richness dependence of Q they found, agrees with our N-body results. Overall, our N-body results of £ are not contradictory 
with all observational results currently obtained. 

(5) These specific but robust features of the cluster three-point correlation function found in this paper, might be an 
important feature of Gaussian gravitational instability theories. Simple toy models of the explosion theory of galaxy formation 
predict a Q which depends rather weakly on v. Therefore a statistical study of £ for redshift samples of rich clusters would 
be very important for testing theories of galaxy formation. 
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